library(dplyr)
library(magrittr)
library(tidyr)
library(ggplot2)

dat <- rio::import("data/Demo_book_country_3.dta")

dem <- dat %>%
  filter(year >= 1789 & year < 2019) %>%
  select(v2x_polyarchy_imp_100, mountains, e_regionpol4, year, country_id)

regs_df <- data.frame(start_year = 1789:1990,
                      end_year = 1818:2019,
                      coef = NA,
                      se = NA,
                      t_value = NA,
                      N = NA,
                      r2 = NA, 
                      series = "lexical")

for (i in 1:nrow(regs_df)){
  subset_data <- dem[which(dem$year %in% c(regs_df$start_year[i]:regs_df$end_year[i])), ]
  reg <- miceadds::lm.cluster(v2x_polyarchy_imp_100 ~ mountains + as.factor(e_regionpol4),
                              data = subset_data, cluster = "country_id")
  reg_sum <- summary(reg)
  regs_df$coef[i] <- reg_sum[2, 1]
  regs_df$se[i] <- reg_sum[2, 2]
  regs_df$t_value[i] <- reg_sum[2, 3]
  regs_df$N[i] <- nrow(reg$lm_res$model)
  regs_df$r2[i] <- summary(reg$lm_res)$r.squared
}

regs_df$year <- regs_df$end_year - 15

g <- regs_df %>%
  filter(series == "lexical") %>%
  ggplot(aes(x = year)) +
  geom_ribbon(aes(ymin = coef - se,
                  ymax = coef + se),
              color = 'lightgray',
              alpha = 0.25) +
  geom_point(aes(y = coef)) +
  scale_x_continuous(breaks = seq(1800, 2025, by = 25)) +
  # ylim(0, 0.04) + 
  labs(x = "Year",
       y = "Coefficient estimates for Mountains") +
  theme(panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.border = element_blank(),
        axis.line = element_line(color = "black"),
        plot.margin = margin(1, 1, 1, 0.25, "cm"))

## ggsave(plot = g, filename = "figure_16_1.png",
##        width = 11, height = 8.5)

ggsave(plot = g, filename = "output/figure_16_1.tiff",
       width = 11, height = 8.5, dpi = 300)
